/* This is a library of useful functions
   that augment those found in math.h
   and the other c libraries.  LRK 06/10 */

#include <stdlib.h>
#include <stdio.h>
#include <math.h>


// returns the minimum value in an array
double getmin(double data[], int N_pts)
{
	int i;
	double min = 10.0e10;
	for(i = 0; i < N_pts; i++) if(data[i] < min) min = data[i];
	return min;
}

// returns the maximum value in an array
double getmax(double data[], int N_pts)
{
	int i;
	double max = -10.0e10;
	for(i = 0; i < N_pts; i++) if(data[i] > max) max = data[i];
	return max;
}

// returns the mean of an array "data"
double getmu(double data[], int N_pts)
{
	int i;
	double sum = 0;
	for(i = 0; i < N_pts; i++)
	{
		sum += data[i];
	}

	return (double)sum/N_pts;
}

// returns the phase a time t for observations
// starting at time t_init with period P
double getphase(double t, double t_init, double P)
{
	double phase;
	phase = (t - t_init)/P - floor((t - t_init)/P);
	return phase; 		
}

//  A histogramming function.  An array "data"
//  is sorted in "hist" by frequency.
void hist(double data[], double hist[], int N_bins, int N_pts)
{
	int  b, i;
	double a;
	
	int min = getmin(data, N_pts);
	int max = getmax(data, N_pts);

	for(i = 0; i < N_bins; i++)
	{
		hist[i] = 0;
	}

	for(i = 0; i < N_pts; i++)
	{
		a = N_bins*(data[i] - min)/(max - min);
		b = (int)a;
		if (b < N_bins && b >= 0) hist[b] = hist[b] + 1;
	}

	for(i = 0; i < N_bins; i++)
	{
		hist[i] = hist[i]/N_pts;
	}	
}

//Stan Eisenstat's getLine() function
char *getLine (FILE *fp)
{
    int c, i;
    int size;                   /* #chars allocated */
    char *line;                 /* Line being read */
    const int DELTA = 16;       /* Initial & incremental storage allocations */

    size = DELTA;
    line = malloc (size * sizeof(*line));
    for (i = 0;  (c = getc(fp)) != EOF && c != '\n';  )  {
	line[i++] = c;
	if (i == size-2) {
	    size += DELTA;
	    line = realloc (line, size * sizeof(*line));
	}
    }
    if (c == EOF && i == 0)  {          /* Check for EOF */
	free (line);
	return NULL;
    }

    if (c != EOF)
	line[i++] = c;
    line[i] = '\0';                     /* Terminate line */

    return (line);
}

//reads one column of double precision data
int readcol(double **a, FILE *infile)
{
	double col;
	int i, n;
	n = 0;
	while(fscanf(infile, "%lf", &col) != EOF) n++;
	*a = (double *)malloc(n*sizeof(double));

	rewind(infile);
	
	i = 0;
	while(fscanf(infile, "%lf", &col) != EOF)
	{
		(*a)[i] = col;
		i++;
	}
	return n;
}

// reads two columns of double precision data
int read2col(double **a, double **b, FILE *infile)
{
	double col1, col2; 
	int i, n;

	n = 0;
	while(fscanf(infile, "%lf", &col1) != EOF && fscanf(infile, "%lf", &col2) != EOF) n++;

	*a = (double *)malloc(n*sizeof(double));
	*b = (double *)malloc(n*sizeof(double));

	rewind(infile);
	
	i = 0;
	while(fscanf(infile, "%lf", &col1) != EOF && fscanf(infile, "%lf", &col2) != EOF)
	{
		(*a)[i] = col1;
		(*b)[i] = col2;
		i++;
	}
	return n;
}

int read3col(double **a, double **b, double **c, FILE *infile)
{
	double col1, col2, col3; 
	int i, n;

	n = 0;
	while(fscanf(infile, "%lf", &col1) != EOF && fscanf(infile, "%lf", &col2) != EOF && fscanf(infile, "%lf", &col3) != EOF) n++;

	*a = (double *)malloc(n*sizeof(double));
	*b = (double *)malloc(n*sizeof(double));
	*c = (double *)malloc(n*sizeof(double));

	rewind(infile);
	
	i = 0;
	while(fscanf(infile, "%lf", &col1) != EOF && fscanf(infile, "%lf", &col2) != EOF && fscanf(infile, "%lf", &col3) != EOF)
	{
		(*a)[i] = col1;
		(*b)[i] = col2;
		(*c)[i] = col3;
		i++;
	}
	return n;
}

int read5col(double **a, double **b, double **c, double **d, double **e, FILE *infile)
{
	double col1, col2, col3, col4, col5; 
	int i, n;

	n = 0;
	while(fscanf(infile, "%lf", &col1) != EOF && fscanf(infile, "%lf", &col2) != EOF && fscanf(infile, "%lf", &col3) != EOF && fscanf(infile, "%lf", &col4) != EOF && fscanf(infile, "%lf", &col5) != EOF) n++;

	*a = (double *)malloc(n*sizeof(double));
	*b = (double *)malloc(n*sizeof(double));
	*c = (double *)malloc(n*sizeof(double));
	*d = (double *)malloc(n*sizeof(double));
	*e = (double *)malloc(n*sizeof(double));

	rewind(infile);
	
	i = 0;
	while(fscanf(infile, "%lf", &col1) != EOF && fscanf(infile, "%lf", &col2) != EOF && fscanf(infile, "%lf", &col3) != EOF && fscanf(infile, "%lf", &col4) != EOF && fscanf(infile, "%lf", &col5) != EOF)
	{
		(*a)[i] = col1;
		(*b)[i] = col2;
		(*c)[i] = col3;
		(*d)[i] = col4;
		(*e)[i] = col5;
		i++;
	}
	return n;
}
	
void bubblesort(double *col1, int n)
{
	int s = 0, t = 0;
	double holder = 0;

	for(s = 0; s < n - 1; s++)
	{
		for(t = s + 1; t < n; t++)
		{
			if(col1[s] > col1[t])
			{
				holder = col1[s];
				col1[s] = col1[t];
				col1[t] = holder;
			}
		}
	}

}
